Tarea 3 Modelo de nicho de la especie Pharomachrus mocinno (Quetzal)

Instalación de paquetes

Esto se hace en la terminal

#install.packages("rJava")
# Paquete para acceder datos en GBIF
#install.packages("rgbif")

# Paquete para acceder datos geoespaciales
#install.packages("geodata")

# Paquete para mapas interactivos
#install.packages("leaflet")

# Paquete para modelado de distribución de especies
#install.packages("dismo")

# install.packages("RColorBrewer")

Carga de paquetes

# Colección de paquetes de Tidyverse
library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.3.3

Warning: package 'ggplot2' was built under R version 4.3.3

Warning: package 'tidyr' was built under R version 4.3.3

Warning: package 'readr' was built under R version 4.3.3

Warning: package 'purrr' was built under R version 4.3.3

Warning: package 'stringr' was built under R version 4.3.3

Warning: package 'forcats' was built under R version 4.3.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# Estilos para ggplot2
library(ggthemes)

Warning: package 'ggthemes' was built under R version 4.3.3

# Paletas de colores de RColorBrewer
library(RColorBrewer)

# Paletas de colores de viridis
library(viridisLite)

Warning: package 'viridisLite' was built under R version 4.3.3

# Gráficos interactivos
library(plotly)

Warning: package 'plotly' was built under R version 4.3.3

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

# Manejo de datos vectoriales
library(sf)

Warning: package 'sf' was built under R version 4.3.3

Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

# Manejo de datos raster
library(terra)

Warning: package 'terra' was built under R version 4.3.3

terra 1.7.83

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

# Manejo de datos raster
library(raster)

Warning: package 'raster' was built under R version 4.3.3

Loading required package: sp

Warning: package 'sp' was built under R version 4.3.3

Attaching package: 'raster'

The following object is masked from 'package:plotly':

    select

The following object is masked from 'package:dplyr':

    select

# Mapas interactivos
library(leaflet)

Warning: package 'leaflet' was built under R version 4.3.3

# Acceso a datos en GBIF
library(rgbif)

Warning: package 'rgbif' was built under R version 4.3.3

# Datos geoespaciales
library(geodata)

Warning: package 'geodata' was built under R version 4.3.3

# Modelado de distribución de especies
library(dismo)

Warning: package 'dismo' was built under R version 4.3.3

# Instalar y cargar el paquete RColorBrewer si no lo tienes
library(RColorBrewer)

Selección de la especie

# Nombre de la especie
especie <- "Pharomachrus mocinno La Llave, 1832"

# Consulta a GBIF
respuesta <- occ_search(
  scientificName = especie, 
  hasCoordinate = TRUE,
  hasGeospatialIssue = FALSE,
  limit = 10000
)

Extracción y guardado de datos

# Extraer datos de presencia
presencia <- respuesta$data

# Guardar los datos de presencia en un archivo CSV
write_csv(presencia, 'presencia.csv')

Lectura de los datos de especie

# Leer los datos de presencia de un archivo CSV
presencia <- read_csv('presencia.csv')

Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
  dat <- vroom(...)
  problems(dat)

Rows: 10000 Columns: 116
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (77): scientificName, issues, datasetKey, publishingOrgKey, installatio...
dbl  (24): key, decimalLatitude, decimalLongitude, crawlId, taxonKey, kingdo...
lgl   (9): isSequenced, isInCluster, eventID, identificationRemarks, organis...
dttm  (6): lastCrawled, lastParsed, dateIdentified, eventDate, modified, las...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

presencia <- st_as_sf(
  presencia,
  coords = c("decimalLongitude", "decimalLatitude"),
  remove = FALSE, # conservar las columnas de las coordenadas
  crs = 4326
)

Gráfico de la distribución por países de la especie Pharomachrus mocinno

# Gráfico ggplot2
grafico_ggplot2 <-
  presencia |>
  st_drop_geometry() |>
  ggplot(aes(x = fct_infreq(countryCode))) +
  geom_bar(
    aes(
      text = paste0(
        "Cantidad de registros de presencia: ", after_stat(count)
      )
    )    
  ) +
  ggtitle("Cantidad de registros de Pharomachrus mocinno por país") +
  xlab("País") +
  ylab("Cantidad de registros de presencia") +
  labs(caption = "Fuente: GBIF") +
  theme_economist()

Warning in geom_bar(aes(text = paste0("Cantidad de registros de presencia: ", :
Ignoring unknown aesthetics: text

# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |> 
  config(locale = 'es')

Gráfico de Distribución porcentual tipo pastel de la especie Pharomachrus mocinno por Provincia en Costa Rica

# Filtrar datos para Costa Rica
presencia_cr <- presencia |> 
  filter(countryCode == "CR")

# Calcular porcentajes de distribución por provincia
distribucion_provincia <- presencia_cr |>
  st_drop_geometry() |>
  group_by(stateProvince) |>
  summarise(count = n()) |>
  mutate(percentage = (count / sum(count)) * 100)

# Filtrar el dataset para eliminar las provincias San José, Guanacaste y Puntarenas
distribucion_provincia_filtrada <- distribucion_provincia |>
  filter(!stateProvince %in% c("Provincia San José", "Provincia Guanacaste", "Provincia Puntarenas"),
         !is.na(stateProvince))

# Crear gráfico tipo pie con ggplot2
grafico_pie <-
  distribucion_provincia_filtrada |>
  ggplot(aes(
    x = "", 
    y = percentage, 
    fill = stateProvince, 
    text = paste0(
      "Provincia: ", stateProvince,
      "<br>Porcentaje: ", round(percentage, 1), "%"
    )
  )) +
  geom_bar(stat = "identity", width = 1, color = "white") +
  coord_polar(theta = "y") +
  ggtitle("Distribución de registros de Pharomachrus mocinno \n por provincia en Costa Rica") +
  geom_text(
    aes(label = paste0(round(percentage), "%")), 
    color = "black",
    size = 3, # Ajustar el tamaño del texto
    position = position_stack(vjust = 0.5) # Para ajustar la posición del texto en cada porción
  ) +
  scale_fill_brewer(palette = "Paired") +  # Cambiar la paleta de colores (puedes probar con otras paletas)
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.title = element_blank()
  )

# Mostrar gráfico
grafico_pie

Mapa de la distribución de la especie Pharomachrus mocinno

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Amanita muscaria"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Registros de Bradypus variegatus"))

Mapa de la distribución de la especie Pharomachrus mocinno junto con variables climáticas de precipitación y temperatura

# Consulta a WorldClim
clima <- worldclim_global(var = 'bio', res = 10, path = tempdir())

# Nombres de las variables climáticas
names(clima)

 [1] "wc2.1_10m_bio_1"  "wc2.1_10m_bio_2"  "wc2.1_10m_bio_3"  "wc2.1_10m_bio_4" 
 [5] "wc2.1_10m_bio_5"  "wc2.1_10m_bio_6"  "wc2.1_10m_bio_7"  "wc2.1_10m_bio_8" 
 [9] "wc2.1_10m_bio_9"  "wc2.1_10m_bio_10" "wc2.1_10m_bio_11" "wc2.1_10m_bio_12"
[13] "wc2.1_10m_bio_13" "wc2.1_10m_bio_14" "wc2.1_10m_bio_15" "wc2.1_10m_bio_16"
[17] "wc2.1_10m_bio_17" "wc2.1_10m_bio_18" "wc2.1_10m_bio_19"

# Definir la extensión del área de estudio
area_estudio <- ext(
  min(presencia$decimalLongitude) - 5, 
  max(presencia$decimalLongitude) + 5,
  min(presencia$decimalLatitude) - 5, 
  max(presencia$decimalLatitude) + 5
)

# Recortar las variables bioclimáticas al área de estudio
clima <- crop(clima, area_estudio)

# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
  # palette = "inferno",
  # palette = "magma",
  palette = rev(brewer.pal(11, "RdYlBu")),
  values(clima$wc2.1_10m_bio_1),
  na.color = "transparent"
)

# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
  # palette = "viridis",
  # palette = "YlGnBu",  
  palette = "Blues",
  values(clima$wc2.1_10m_bio_12),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c("Temperatura", "Precipitación", "Registros de Pharomachrus mocinno")
  ) |>
  hideGroup("Precipitación")

Creación del modelo de nichos para la especie Pharomachrus mocinno

# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
  decimalLongitude = presencia$decimalLongitude,
  decimalLatitude = presencia$decimalLatitude
)

# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)

# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)

# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)

# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7) 
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
  1:n_presencia, 
  size = round(0.7 * n_presencia)
)

# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]

# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]

# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a es este el que acepta el paquete dismo
clima <- raster::stack(clima)
# Ejecutar el modelo
modelo_maxent <- maxent(x = clima, p = entrenamiento)

Loading required namespace: rJava

# Aplicar el modelo entrenado a las variables climáticas 
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, clima)

# terra::extract() extrae los valores del raster de predicción 
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos 
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
  prediccion, 
  evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)

# Generar puntos aleatorios dentro del área de estudio definida. 
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = clima, n = 1000)

# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
  prediccion, 
  ausencias
)

# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)

Curva ROC y resultado del dato AUC del modelo par ala especie Pharomachrus mocinno

# Datos para graficar la curva ROC
datos_roc <- data.frame(
  FPR = resultado_evaluacion@FPR,
  TPR = resultado_evaluacion@TPR,
  Umbral = resultado_evaluacion@t
)

# Valor AUC
auc <- resultado_evaluacion@auc

# Gráfico ggplot2
grafico_ggplot2 <-
  ggplot(
    datos_roc, 
    aes(
      x = FPR, 
      y = TPR,
      u = Umbral
    )
  ) +
  geom_line(
    color = "blue", 
    size = 1
  ) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
       x = "Tasa de falsos positivos (FPR)",
       y = "Tasa de verdaderos positivos (TPR)") +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

# Gráfico plotly
ggplotly(grafico_ggplot2) |> 
  config(locale = 'es')

Mapa interactivo de idoneidad para la probabilidad de encontrar la especie Pharomachrus mocinno

# Paleta de colores del modelo
colores_modelo <- colorNumeric(
  palette = c("white", "black"),
  values(prediccion),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>  
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage( # capa raster de temperatura
    clima$wc2.1_10m_bio_1,
    colors = colores_temperatura, # paleta de colores
    opacity = 0.6,
    group = "Temperatura",
  ) |>
  addRasterImage( # capa raster de precipitación
    clima$wc2.1_10m_bio_12,
    colors = colores_precipitacion, # paleta de colores
    opacity = 0.6,
    group = "Precipitación",
  ) |>
  addRasterImage( # capa raster del modelo de distribución
    prediccion,
    colors = colores_modelo,
    opacity = 0.6,
    group = "Modelo de distribución",
  ) |>  
  addCircleMarkers(
    # capa de registros de presencia (puntos)
    data = presencia,
    stroke = F,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>  
  addLegend(
    title = "Temperatura",
    values = values(clima$wc2.1_10m_bio_1),
    pal = colores_temperatura,
    position = "bottomleft",
    group = "Temperatura"
  ) |>
  addLegend(
    title = "Precipitación",
    values = values(clima$wc2.1_10m_bio_12),
    pal = colores_precipitacion,
    position = "bottomleft",
    group = "Precipitación"
  ) |>
  addLegend(
    title = "Modelo de distribución",
    values = values(prediccion),
    pal = colores_modelo,
    position = "bottomright",
    group = "Modelo de distribución"
  ) |>  
  addLayersControl(
    # control de capas
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Temperatura",
      "Precipitación",
      "Modelo de distribución",
      "Registros de Pharomachrus mocinno"
    )
  ) |>
  hideGroup("Temperatura") |>
  hideGroup("Precipitación")

Mapa de idoneidad binario para la especie Pharomachrus mocinno

# Definir el umbral
umbral <- 0.5

# Crear el raster binario
prediccion_binaria <- (prediccion >= umbral) * 1

# Crear la paleta de colores para el raster binario
colores_prediccion_binaria <- colorFactor(
  palette = c("transparent", "blue"),  # "transparent" para las áreas no adecuadas
  domain = c(0, 1),
  na.color = "transparent"
)

# Mapa
leaflet() |>
  addTiles(group = "Mapa general") |>
  addProviderTiles(
    providers$Esri.WorldImagery, 
    group = "Imágenes satelitales"
  ) |>
  addProviderTiles(
    providers$CartoDB.Positron, 
    group = "Mapa blanco"
  ) |>
  addRasterImage(
    prediccion_binaria,
    colors = colores_prediccion_binaria,
    opacity = 0.6,
    group = "Modelo de distribución binario",
  ) |>
  addCircleMarkers(
    data = presencia,
    stroke = FALSE,
    radius = 3,
    fillColor = 'red',
    fillOpacity = 1,
    popup = paste(
      paste0("<strong>País: </strong>", presencia$country),
      paste0("<strong>Localidad: </strong>", presencia$locality),
      paste0("<strong>Fecha: </strong>", presencia$eventDate),
      paste0("<strong>Fuente: </strong>", presencia$institutionCode),
      paste0("<a href='", presencia$occurrenceID, "'>Más información</a>"),
      sep = '<br/>'
    ),
    group = "Registros de Bradypus variegatus"
  ) |>
  addLegend(
    title = "Modelo de distribución binario",
    labels = c("Ausencia", "Presencia"),
    colors = c("transparent", "blue"),
    position = "bottomright",
    group = "Modelo de distribución binario"
  ) |>
  addLayersControl(
    baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
    overlayGroups = c(
      "Modelo de distribución binario",
      "Registros de Pharomachrus mocinno"
    )
  )

Warning in colors(.): Some values were outside the color scale and will be
treated as NA